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Alan D. Freed 
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Abstract 

Higher-level physical laws applicable to biological tissues are presented that will 
permit the modeling of metabolic activity at the cellular level, including varia- 
tions in the mass of a tissue. Here the tissue is represented as a fluid/solid mixture, 
wherein molecular solutes transport within the fluid, and cells can migrate through- 
out the porous solid. Variations in mass can arise via exchanges in mass between the 
constituent phases within a control volume such that mass is conserved in the tis- 
sue overall. The governing balance laws for mass, momentum, energy, and entropy 
are a special case of those describing a chemically reacting mixture with diffusion. 
Thermodynamic constraints on the constitutive structure are addressed. 


1. Introduction 

In the absence of gravity, or in reduced gravity environments like on the surfaces 
of the Moon or Mars, astronauts’ bodies will undergo physical changes brought 
about by these varied states of gravity (Buckey 2006). The forces carried by their 
musculoskeletal frames will be altered, because their weights will be reduced. Many 
vital tissues will resorb and remodel due of these changes in gravity, altering their 
masses, densities, and architectures. This will bring about changes in the physi- 
cal capabilities of these astronauts, with the potential of placing their health at 
increased risk (compared to their health risk here on Earth). If an injury were to 
occur, it could increase the health risk of the remaining crew, and even impact 
the ability of the crew to achieve its mission objectives. Higher-level physical laws 
pertinent to investigating such risk scenarios are presented in this paper. 

To address many of the specific risk scenarios pertaining to astronaut health 
during long duration voyages deep into space will require numeric simulations, as 
no data base from past experiences exists to draw inference from at this time. These 
simulations will likely be based upon biochemical, biomechanical, and mechanobio- 
logical models that span the length and physiologic scales of molecules, cells, tissues, 
organs, and organisms. We are of the opinion that studies intended to bridge such 
vast length and physiologic scales will require the cells, solute molecules, interstitial 
fluid, and extra cellular matrix of living tissue to be addressed individually in an 
unifying framework. The physical laws presented in this paper provide a theoretical 
foundation for studies such as these to be done. 

Relevant fields of study include (G. A. Holzapfel 2006, personal communication): 
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• Biophysics: The science of biology and medicine, as revealed through appli- 
cations of physical laws and theories. 

• Biomechanics: The study of living systems through developments, extensions, 
and applications of mechanics targeting to better explain phenomena in biol- 
ogy, medicine, and bioengineering. Biomechanics focuses on whether or how 
function follows structure via physical laws. 

• Mechanobiology: The study of biological reactions of cells in response to 
changes in their mechanical environment, as in growth, remodeling, adapta- 
tion, and repair. Mechanobiology focuses on whether or how structure follows 
function via biological laws. 

This is a biophysics paper whose future applications will reside primarily in the 
discipline of mechanobiology and, to a lesser degree, in biomechanics. 

In an assessment of the role of computational sciences in the development and 
application of higher-level physical laws, like mixture theory, the Nobel Prize Lau- 
reate for Physics in 1998, Robert Laughlin, recently stated the truism (Laughlin 
2002): “It is not generally possible to start from the wrong equations and get the 
right result.” Our intention is to pay particular attention in getting the physics 
correct, within the confines of a predefined set of assumptions that pertain to bi- 
ology, so that numeric computations conforming with these assumptions will have 
the capability of making accurate predictions. 


(a) Mixture Theories in Tissue Mechanics 

Binary models of porous elastic solids (cells, collagen, elastin, bone, muscle, 
etc.) saturated with inviscid fluids (water, saline, plasma, lymphatic fluids, etc.) go 
by the names of ‘biphasic theory’ in the soft-tissue mechanics literature (cf. Mow 
et al. 1980) and ‘poroelasticity’ in the hard-tissue mechanics literature (cf. Cowin 
2001). Both are mixture theories that can trace their origins back to Biot’s (1941) 
original theory for soil consolidation, to which Truesdell (1957) gave rigor. They are 
equivalent theories under certain conditions, most notably, under an assumption of 
intrinsic incompressibility (Kenyon 1978). 

A ‘triphasic theory’ was developed latter on to improve upon the predictions 
of biphasic theory in the osmotic swelling of cartilage. This was achieved through 
the actions of charged anions and cations introduced as solute phases into their 
mixture formulation (Lai et al. 1991). ‘Bicomponent theory’ is a simpler alternative 
to multi fluidic-phase theories in that it has treated filtration, osmotic swelling, 
and buoyancy effects as separate local forces of interaction between the liquid and 
solid components, without the need to introduce additional fluidic phases (Lanir 
1987). A recent comparison between the triphasic and bicomponent theories has 
demonstrated that they possess like capabilities (Wilson et al. 2005). 

Another modification to mixture theory was introduced by Mak (1986), who 
replaced the elastic matrix of biphasic theory with a viscoelastic matrix in order to 
account for the viscoelastic attributes of the intrinsic gels (proteoglycans, hyaluronic 
acids, etc.) and matrix constituents (collagen, muscle, etc.) that convect with the 
porous solid in an assumed affine manner. 
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(b) Growth Theories in Tissue Mechanics 


Cells will change the mass, architecture, and/or volume of both themselves 
and their surrounding tissue, the extra cellular matrix, in their effort to maintain 
homeostasis (Wang & Thampatty 2006). Growth, in a general usage of the word, 
is said to occur whenever the volume fraction of the porous solid increases at the 
expense of the volume fraction of the interstitial fluid due to cell activity. 

There has been a recent flurry of papers addressing the physical laws of mass, 
motion, and thermodynamics pertaining to open-system continua where mass is 
being created locally (e.g., Epstein & Maugin 2000; Klisch & Lotz 2000; Humphrey 
& Rajagopal 2002; Lubarda & Hoger 2002; Klisch & Hoger 2003; Kuhl & Steinmann 
2003; Garikipati et al. 2004; Menzel 2005; and Guillou <fc Ogden 2006). These studies 
extend the earlier works of Hsu (1968), Cowin & Hegedus (1976), and Skalak et al. 
(1982). A common class of problems mentioned in all of these papers is tissue 
growth. 

Cowin k, Hegedus (1976) derived balance laws for a porous single-phase contin- 
uum. Klisch & Lotz (2000) applied balance laws for a fluid/ solid mixture. Humphrey 
& Rajagopal (2002) applied balance laws for a fluid/multiple-solid mixture, where 
they constrained all solid constituents to deform with the same affine motion, and 
then applied the rule-of-mixtures from homogenization theory to quantify stress in 
the solid. Klisch & Hoger (2003) provide balance laws for tissues derived from many 
mixture scenarios. Garikipati et al. (2004) derived balance laws for a multiple-fluid/ 
solid mixture. And the others applied balance laws that pertain to a dense single- 
phase continuum. Our approach is to consider a four-phase mixture that we show to 
be a special case of Truesdell’s (1957) theory for chemically reacting mixtures with 
diffusion. Here, the fluid transports a population of solutes needed for metabolism, 
and the solid hosts a population of cells that utilize these metabolites. 

All of the above mentioned tissue theories allow for volumetric growth, while the 
theories of Epstein & Maugin (2000), Kuhl & Steinmann (2003), Garikipati et al. 
(2004), and Guillou & Ogden (2006) allow for an additional growth through a flux 
of mass. In contrast, our theory introduces growth and other biological activities 
through various mass exchanges between the constituents present within a mass 
element. 

Most of these theories introduce a multiplicative decomposition of the deforma- 
tion gradient, e.g., F = F e Fg (cf. Taber 1995), wherein the elastic F e and growth 
F g tensors are often interpreted with slight variations betwixt them, which usually 
include some form of further multiplicative refinement. In such theories, growth 
is a tensor field, viz. Fg, whose evolution must be specified by some constitutive 
law. Such a decomposition of the deformation gradient is always permissible in 
single-phase continua. We do not introduce a multiplicative decomposition of the 
deformation gradient; we cannot, because F is not defined in our theory. However, 
one could define such a decomposition for the extra cellular matrix in our mixture 
theory, viz., F m = F"'F™, if one so desired. 

An alternative approach to the notion of a growth tensor is found in Humphrey 
& Rajagopal (2002), where they introduced the idea of a survival function as a 
means for assigning natural configurations to constituents belonging to mass points 
at the time of their creation. This function accounts for the rate of production and 
the life span of each solid constituent. It enters into the theory as a kernel function 
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in a Volterra integral equation of the second kind, like those used in the literature of 
viscoelastic liquids. Of such configurations, they state: “It is differences in natural 
configurations, not incompatible strains, that likely leave a tissue stressed in the 
absence of applied tractions.’’ We find this to be a more appealing approach than 
that of decomposing the deformation gradient, as it better reflects the interface 
between physics and biology. 


(c) Approach Taken 

A volume element of tissue is considered to be comprised of four types of con- 
stituents: an interstitial miscible liquid (iml) suspending a population of solutes 
(pos), and an extra cellular matrix (ecm) hosting a population of cells (poc). Each 
constituent may itself be an assembly of many other sub-constituents, but that level 
of sophistication is not required at this time in order for us to be able to derive a set 
of higher-level physical laws that govern tissue and permit the modeling of biolog- 
ical activity. However, such distinctions will become necessary when constructing 
mechanobiologic constitutive models for specific tissues, which lies beyond the scope 
of the present paper. 

The various documents mentioned above, where physical laws have been de- 
rived for tissue growth, are applicable to open-system continua; that is, mass can 
be created and/or destroyed within a material element in these theories. The the- 
ory presented here describes a closed-system continuum where mass is conserved 
within a material element. It is a fluid/solid mixture theory, where the fluid is com- 
prised of an iml suspending a pos, while the solid is comprised of an ecm hosting 
a poc. Cells manage themselves and their surrounding environment by, in essence, 
converting mass from one form into another in order to suit some function, be it 
biological, chemical, or physical (Wang & Thampatty 2006). As an analog, cells are 
factories that, in a physical sense, regulate exchanges in mass between the various 
constituents of tissue. The nutrients that these cellular factories run on arrive as 
molecules suspended by the liquid that surrounds them; furthermore, their wastes 
depart as suspended particles in this same bathing liquid. So, in our approach to 
mechanobiology, mass is neither created nor destroyed at a continuum point, but 
rather, it is moved from one phase to another through biological processes. 

Adopting the premise of Lai et al. (1991), the fluid phase of tissue is considered to 
be comprised of a miscible liquid (the solvent or iml) that suspends soluble molecules 
(the solutes or pos). In like manner, adopting the simplifying premise of Humphrey 
& Rajagopal (2002), the solid phase (the porous lattice or ecm) is considered to be 
comprised of various constituents that experience the same affine deformation, each 
with their own natural configuration. Managing this structure is the responsibility 
of a poc. If the stresses carried by the various solid constituents are assumed to 
sum according to the rule-of-mixtures, as dictated by homogenization theory, then 
tractions can be assigned at the boundaries in a straightforward manner. These 
assumptions are idealizations, and exceptions can be put forward for each of them. 
Nevertheless, many important problems can be solved where these assumptions 
hold, and it is the physical laws that govern this class of problems that we choose 
to address. 

The scope of our theory is quite broad in that it will allow for the modeling of 
cellular activity (e.g., cell birth, death, metabolism, and migration), tissue behavior 
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(e.g., growth, swelling, remodeling, healing, and nutrient /waste transport), and 
mechanical response (e.g., stress, strain, and energy), all within the same framework. 
Applications of our theory to specific tissues undergoing specified biological and/or 
physical processes are left to future papers. 

Supposition 1 . Biological tissues are porous materials saturated with miscible fluids 
that suspend solute molecules. The solid material is comprised of an assembly of 
constituents that experience the same affine deformation, but that possess different 
natural configurations. This porous solid contains a distribution of cells that attach 
to the matrix and are capable of locomotion by migrating through its cavities. Cells 
are the caretakers of tissue, and in the execution of their various duties, they con- 
vert mass from one form belonging to some constituent into another form possibly 
belonging to a different constituent in a manner that conserves mass overall. 
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2. Terminology 

(a) Acronyms 


extra cellular matrix a = m poc population of cells a = c 

interstitial miscible liquid ..a=l P os population of solutes a = s 


(b) Nomenclature 


free-energy supplies, a e 
{c,l,m,s} 
volume in 3-space 
density of mass supplies 
material derivatives 
differential element of area 
differential elements of mass 
differential elements of volume 
diffusion coefficient of a cell type 
solute driven cell diffusion coeffi- 
cient 

diffusion coefficient of a solute 
type 

strain-rate tensors 

density of internal-energy supplies 

basis vector, i = 1,2,3 

deformation gradient tensors 

gravity vector 

time-step size 

identity tensor 

Jacobian of the ecm deformation 
Boltzmann’s constant 
velocity gradient tensors 
molecular weight of solute j 
molar concentration of solute j 
outward unit normal vector 
number of cell species 
number of solute species 
momentum supplies 
heat flux vectors 


r 

radius 

r, r a 

heat productions 

d“ 

density of entropy supplies 

t 

time 

t, t a 

traction vectors 

T 

Cauchy stress tensor 

r pO! 

partial stress tensors 

u a 

diffusion velocity vectors 

v a 

velocity vectors 

W, W“ 

vorticity tensors 

X 

current position vector 

X 

velocity vector 

X 

acceleration vector 

X m 

reference vector for the ecm 


arbitrary mass-averaged fields 

y, y« 

density of entropy productions 

8£ 

surface in 3-space 

e, e a 

internal energies 

r> 

viscosity- 

rh rf 

entropies 

9 

absolute temperature 

li* 

mass fractions 

P, P a 

partial mass densities 

P% 

mass transfer rate from to a 

Q a 

actual mass densities 

4 >“ 

volume fractions 

<P 

motion map for the ecm 

4/, \jf a 

Helmholtz free energies 

i/r 

motion map for the iml 
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3. Mixtures 


We adopt the approach of mixture theory (Truesdell 1957) to implement our sup- 
position, and to establish the physical laws that govern it. The reader is referred 
to the review article by Bowen (1976), the text by Rajagopal & Tao (1995), and 
the references cited therein for more thorough treatments on mixture theory. The 
former gives an excellent presentation of the classical theory for mixtures, while the 
latter presents a nice collection of solved boundary-value problems. 

The derivation of constitutive equations for specific tissue types, and the con- 
struction of variational principles required to solve boundary value problems that 
pertain to these constitutive formulae are the topics of future research endeavors, 
and their publications. 


(a) Masses and Volumes 

At any given instant in time, let a differential element of volume be comprised 
of four distinct sub- volumes 

da = y^da“ = da c + d?/ + Av m + da s , a e {c,£,m, s}, (3.1) 

where properties with an index of ^ or s designate an affiliation with the iml 
or pos phases of the fluid, and properties with an index of m or c designate an 
affiliation with the ecm or poc phases of the solid, respectively, while an index of “ 
can represent any one of these four constituents, and a summand without limits is 
considered to sum a over all four of its phases. Volume fractions are then defined 
as 

(p a — da“/da with 0 < <p a < 1 = 1. (3.2) 

which sum to 1 because tissues are considered to be simply connected. The volume 
fraction of each constituent will vary between tissue type. In the literature, the 
volume fraction of fluid, which in our theory is quantified by <j + 0 s , is called the 
porosity. 

The total mass dm of a volume element da is the sum of its constituent masses 
in that 

d /77 = d m a . (3.3) 

The density of this mass element is 

p = din/dv, (3.4) 

which has partial mass densities p a that satisfy 

p — wherein p a = dm a /dv, (3.5) 

and actual mass densities g a that are given by 

Q a = dm a /dv a with p a = (\> a Q a p = J2<P a e c ‘- (3-6) 

It is also useful to define mass fractions (or concentrations) as 

p, 01 — dm “/dm = p a /p y^ /i“ = 1, (3.7) 

which sum to 1 because tissues are saturated mixtures. 
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(b) Kinematics 

We select a Cartesian coordinate system (1,2, 3) described by a set of orthonor- 
mal base vectors {ei,c 2 , ^3 }, where a ■ ej = 8ij and a x ej = eijk^k, with the 
repeated index k being summed from 1 to 3 in the usual manner. Here • and x 
denote the inner and cross products, respectively, while Sy signifies the Kroneker 
delta, and eijk is the permutation symbol. 

At current time t, let all four constituent mass points dm“ co-habitate with one 
another at some volume element di> resulting in a mass element dm located by the 
coordinates {x\,X 2 ,x-}}, or equivalently, by the position vector x = Xj ej. 

(i) Primitive Variables 

The velocity vector is typically selected as the primitive kinematic variable in 
fluid mechanics; whereas, in mixture theory, position vectors belonging to the vari- 
ous phases are usually assigned as its primitives. Here we break with this tradition 
and select velocity vectors as the primitive kinematic variables for the iml, poc & 
pos, and we employ a position vector as the primitive kinematic variable for the 
ecm. This mitigates the need to establish reference configurations for all constituents 
except for the ecm. 

In accordance with supposition 1 , we assign a set of coordinates {X™, X™, X™ } to 
a mass point dm m of an ecm to establish its reference location X m — X™ at some 
initial time to- Let the motion of this mass point through space be described by 

X m = (p~ 1 (x,t) and x = (p{X m ,t), (3.8) 

where the vector mapping function (p = cp/ ej is considered to be continuous, suffi- 
ciently differentiable, and invertible; hence, (p~ l exists. The velocity of an ecm mass 
point is therefore given by 

v m (x,t) = dfcptet, (3.9) 

wherein 9 f (») = 9(»)/9? is the partial derivative taken with respect to time t, with 

• denoting an arbitrary field. 

Also in accordance with supposition 1 , let the motion of mass point dm 1 for the 
iml be described by 

v i = Tjr{x,t), (3.10) 

where the vector mapping function ifr = quantifies the velocity of the miscible 
liquid, and is considered to be continuous and sufficiently differentiable, but unlike 
(p, jr need not be invertible. 

Collectively, the iml & ecm constitute a classic fluid/solid mixture to which we 
now add the attributes needed to incorporate biology. 

In accordance with supposition 1, let the motion of mass point dm s for the pos 
be described by the vector mapping 

N s N s 

v s = f (x,t) — ^ £>/ grad In p] with p s = ^ pj , (3.11) 

i = 1 i = 1 

where &■ is the solute diffusion coefficient, and p. is the mass density, both be- 
longing to species i in some pos that contains N s species. The gradient operator 
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grad(*) = 3(«)/3 jc,- £,■ denotes a partial derivative taken with respect to the current 
position x. The first term on the right-hand side of equation (3.11) describes the 
fluid velocity of the liquid that suspends these particles as it flows through the 
tissue; it is the velocity map of the iml given in equation (3.10). The second term 
accounts for a perturbation to this mean velocity field caused by Brownian mo- 
tion of the solute particles, as described by Fick’s first law for random diffusion. 
The minus sign ensures that solutes migrate from regions of high concentration (or 
density) to regions of low concentration. 

Einstein’s theory for Brownian motion, the topic of his 1905 Ph.D. thesis, pro- 
vides a formula for quantifying the diffusion coefficient of a spherical body with 
radius r suspended in a fluid of viscosity r] that is subject to the random bombard- 
ment of atoms from within the fluid; it being, 

5) = k6/6nrjr , (3.12) 

where k is Boltzmann’s constant, and 9 is the absolute temperature (i.e., 9 > 0). 
This formula can be applied to quantify, for example, the diffusion coefficients of 
globular proteins, to which the £)■ refer. However, this formula must not be used 
for purposes of quantifying the diffusion coefficients associated with cell migration 
presented below. These coefficients need to be experimentally determined. The pos 
diffuse according to Einstein’s physics for Brownian motion; whereas, the poc diffuse 
under their own locomotive processes. The diffusion coefficient 3) s will diminish 
whenever the solute size approaches the pore size of the matrix material, which 
varies with tissue type, and will become zero whenever the solute size exceeds the 
pore size, with the tissue now acting as a filter. 

Finally, in accordance with supposition 1, let the motion of mass point d m c for 
the poc be described by the vector mapping 

N c N c N s N c 

v c = d t (p~Y J £>i grad In Pi +J2H Dfj grad In p S j with / = £*>?• (313) 

i — 1 i=l 7 = 1 1 = 1 

where £)f is the cell diffusion coefficient associated with random cell locomotion, 
and p\ is the mass density, both of species i belonging to some poc that contains N c 
species, while 3) ff is the cell diffusion coefficient of cell species i whose locomotion 
occurs as a response to the presence of solute species j . The first term on the right- 
hand side of equation (3.13) is a passive contribution that describes the velocity of 
the deforming matrix to which the cells adhere; it is the time derivative of the ecm 
motion map given in equation (3.8). The second term is an active contribution that 
accounts for the haptotaxis migration of cells along a cell concentration gradient, 
which is described by Fick’s first law for random diffusion in accordance with, e.g., 
Nobel’s (1987) experimental observation that cell locomotion is Markovian. The 
third term is another type of active contribution to cell diffusion. This term accounts 
for the chemotaxis migration of cells along a chemical concentration gradient. The 
positive sign affiliated with this term implies that the diffusion coefficient will be 
positive whenever cells are attracted to a chemical stimulus, and negative whenever 
they repel from it. Not all cell species are capable of locomotion, so some of the 3)f 
will be zero valued, and as such, their S)ff will be zero valued, too. Furthermore, 

‘■J 

most solute species j in a pos do not serve as an attraction potential for cells of 
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species i that are capable of locomotion and which belong to some poc; hence, most 
of the coefficients S)ff will be zero valued for cell species i where S)f is non-zero. 

The physical laws derived herein, as they apply to the pos and poc, pertain to 
their collective populations, not to the individual species of which they are com- 
prised. Equations (3.11 & 3.13) quantify averaging schemes used to establish these 
kinematic responses. 


(ii) Kinematic Variables 

Mixture theory averages the velocities of its constituent phases, which we sup- 
pose obey equations (3.9-3.11 & 3.13), in order to obtain a representative velocity 
field for the mixture which, in our case, we conjecture to be descriptive of tissue. 
In particular, the mean (or barycentric) velocity is defined as 

x = ^/z“ v a or equivalently px — ^ ^p a v a , (3.14) 

which is a mass- weighted average of the local velocity fields (Truesdell 1957). The 
diffusion velocity 

u a — v a — x .-. ^p“h“ = 0, (3.15) 

proves to be a useful measure describing the velocities of the constituents relative 
to the overall velocity of the mixture. 

By defining material derivatives for the phase constituents as 

d“ (•) = 9 f (») + grad (•) • v a , (3.16) 

where d“ (•) denotes the total derivative taken with respect to time t of a field 
belonging to the a th phase, one obtains a material derivative for the barycentric 
frame that is given by 

d?( # ) = 9 f (») + grad (•) • i pd,(«) = ^p“d“(«), (3.17) 

where d f (») denotes the total derivative taken with respect to time t of a field 
belonging to the continuum mixture, such that 


d“ (•) = d f f •) + grad (•) • u a , 


(3.18) 


which is another identity of value. It is vital to distinguish between the material 
derivatives of the phases d“ from the material derivative of the mixture d f during 
the construction of a mixture theory, and in its applications. 

The ecm has a deformation gradient tensor associated with it that is defined by 


F m = GRAD <p 
(F m ) -1 = grad<p -1 


d(pi 

dx T ‘ 

d<Pi l 


dxj 




e J' 


(3.19) 


so that F m (F m ) 1 = (F m ) 1 F m — I, wherein I = Sy ei <8> e/ is the identity tensor, 
with operator <g> denoting the dyadic product, which satisfies (e/<8> ej) v = (v • ej) et 
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for all vectors v. Because of the invertibility of mapping <p, it follows that the 
Jacobian of this deformation is non-zero, viz., 

J m = detF m /0, (3.20) 

where det designates a matrix determinant, in this case, of components F = 
d(pi/dX'j l where F m = F"j ei <g> ej. Because velocity is the primitive variable in 
the other three phases, these phases do not associate with a reference configuration 
a priori, and as such, they do not possess deformation gradient tensors that can be 
constructed in said manner. We will return to this point later. 

Velocity gradient tensors for the four separate phases are defined by 

O Oi 

L“ = grad v a — ej® (3.21) 

axj 

whose associated velocity vectors are given in equations (3.9-3.11 & 3.13), and 
whose components are L “. = 3v“/3 Xj. The mean velocity gradient is similarly de- 
fined by 

3 Xi 

L = grad x = - — ej® e ; , (3.22) 

OXj 

which, incidently, relates to the velocity gradients of the constituents through the 
formula 

pL — ^(p“L“ + u a ® gradp“). (3.23) 

a 

Hence, equations (3.14 & 3.23) establish maps between kinematic fields of the sub- 
continua (i.e. , the phases) and the continuum (i.e., the mixture). 

The symmetric parts of the constituent and mean velocity gradients, 

D“ = i(L“ + (L“) t ) and D = I(L + L T ), (3.24) 

define their strain rates, while their associated skew-symmetric parts, 

W“ = £(L“ - (L°') T ) and W = i(L -L T ), (3.25) 

quantify their vorticities, and as such 

L“ = D“ + W" and L = D + W, (3.26) 

wherein T designates the transpose, e.g., given L = Uj ej®ej then L T = L yi - ej®ej. 

(iii) Deformation Gradients from Velocity Gradients 

The local kinematic fields relate to the global kinematic fields through the veloc- 
ities (equation 3.14) and velocity gradients (equation 3.23). In stark contrast, local/ 
global mappings do not exist for either an initial position vector X or a deformation 
gradient tensor F, regardless of whether one chooses velocities or displacements as 
the kinematic primitives (cf. Bowen 1976). 

Simple fluids with memory (Coleman & Noll 1964), for example, require a knowl- 
edge of the deformation history of the fluid; in particular, of the relative deformation 
gradient F f (r) = 3 Xj(t)/dxj(r) e ,• <8> ej defined for any reference time r, r < t. In 


NASA/TM— 2007-214827 


10 



tissue mechanics, relative deformation gradients of a solid phase have been used by 
Humphrey & Rajagopal (2002) to handle the birth and death events of the con- 
stituents therein. To obtain such measures of deformation for the poc, iml & pos 
phases in our formulation, where the velocity vectors v c , v e & V s are the primitives, 
not the position vectors x c , x i & X s , will require an integration of their gradients. 

One can begin by discretizing the interval of integration [0, t ] into N sub-intervals 
such that 0 — to < t\ < ■ ■ ■ < t^-i < tjy = t. Provided that the velocity gradient 
L“ is known, and its associated deformation gradient F“ is being sought, as in a 
finite element implementation of our theory, one can start by assigning an initial 


condition of 

F,“(**) = I, (3.27) 

so that F“ (tk) can be approximated for any n > k, given that k = 0, 1, 2 n — 1 

and n — 1,2,..., iV, by employing, e.g., a mid-point predictor 

F“ (t k ) = F l_ 2 (t k ) + 2AL“(f JI _ 1 )F“_ 1 (f*), n > k + 2, (3.28) 

that can be started at n — k + 1 with the forward-Euler predictor 

F^ +1 fe) = I + ^L“fe), (3.29) 

which are then corrected with the trapezoidal rule 

F lit k ) = F “_ 1 (f*) + (3.30) 

when advancing the solution over an uniform time step of h = t n — t„-\ V «. Here 
we are solving the governing differential equation F“ = L"F“ using a predictor/ 
corrector with accuracy 0(/? 3 F“). Other integrators could be used, too. 

4. Mass Balance in Tissues 

Any good textbook in continuum mechanics can be consulted to acquire detailed 


derivations of the conservation laws that govern continua (e.g., Holzapfel 2000). 
Here we use axioms for stating the physical laws that govern continua, and pos- 
tulates for stating extensions to these axioms that are assumed to apply to the 
sub-continua, viz., the phases of a mixture. 

Consider a connected and bounded region fixed in 3-space in an Eulerian 
frame at current time t that is enclosed by a surface 8£. Such a region establishes 
a control volume through which a physical law can be transcribed from axiom into 
formula. 


(a) Balance of Mass 

Axiom 1. The rate at which mass increases inside of 3 equals the flux of mass 
entering across 83. 

Axiom 1 is a statement for mass conservation in a continuum, which assumes 
the form of an integral equation; in particular, 

3, / pdv — — / px-n&a, (4.1) 

Js JstB 
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where dt> is an element of volume, da is an element of surface area, and n is its 
unit normal. The surface integral is negated due to the convention that n points 
outward. Provided that the integrands are continuous and sufficiently differentiable, 
then an application of the divergence theorem to the surface integral in equation 
(4.1) allows the resulting integral equation to be recast as the field equation 

d t p + pdivi = 0, (4.2) 

with div v = dvi/dxi signifying the divergence operator, and where use has been 
made of the identity 

div(pv) = grad(p) • v + pdivto (4.3) 

Our theory differs from the existing theories for tissue growth discussed in the 
Introduction of this paper in that here mass is conserved. Our theory describes 
a closed-system continuum; whereas, prior growth theories describe open-system 
continua, viz., their counterparts to our equation (4.2) have values other than zero 
on their right-hand sides. 

Postulate 1. For a phase in a mixture, the rate at which mass increases inside 
of 3 equals the flux of mass entering across 83, plus any masses that are being 
exchanged between it and the other phases within 3. 

As a matter of notation, let /6^*)”, atl0n represent the rate at which mass is being 
moved from a ‘source’ phase to a ‘destination’ phase within a control volume. With 
there being four separate phases, there exist twelve possible exchange rates. Only 
a few are expected to be active for any given tissue type; nevertheless, all twelve 
cases are considered in the construction of our general theory. What functional 
forms these terms may take on is tissue dependent. This is a constitutive modeling 
issue, and as such, is not addressed herein. It is through these functions, in part, 
that biological laws and physical laws can interact. 

The mass balance governing the ecm within a tissue, as described in supposi- 
tion 1 and constrained by postulate 1, obeys the integral equation 

9 1 f p m dv = - f p m v m ■ n d a + I (p m c - p, c n + p m t - + p m s - p s m ) da, (4.4) 

Jb Jsb Jb 

whose local form is the field equation 

d™ p m + p'"div v m — c m with c m = P m c ~ P c m + P m t - pl + P™ ~ P s m , (4.5) 

where c m is the density of mass supply to the extra cellular matrix, i.e., the rate 
per unit volume at which mass is being moved from the three other phases into the 
ecm, minus the rate at which it is being removed from the ecm and placed into the 
other three phases. This formula can be rewritten as 

d?(p m detF m ) = « w detF m , (4.6) 

because of the identity d^detF'" = detF'"divu'" (cf. equation 1.2.12 in Bowen 
1976), from which it immediately follows that 

detF m = 1. (4.7) 
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whenever the matrix phase is incompressible, which is a good assumption for soft 
tissues, but not for hard tissues. 

The mass balance governing the poc within a tissue, as described in supposition 1 
and constrained by postulate 1, obeys the integral equation 


3 , 


f P c dv = - f , 
Ji 8 J&B 


p c v c ■ nda + 


L {i, ‘ 


Pi + Pm 


+ Ps 


Pc) dth (4.8) 


whose local form is 


d c t p c + p c div v c = € c with = f>\ - p[ + p c m - p m c + p c s - p s c , (4.9) 

where c c is the density of mass supply to the population of cells, i.e., the rate per 
unit volume at which mass is being moved from the three other phases into the 
poc, minus the rate at which it is being removed from the poc and placed into the 
other three phases. From this equation, it immediately follows that 

trL c = div v c = 0, (4.10) 


whenever the cells are incompressible, which is a good assumption, with trL = L a 
signifying the trace operator. 

The mass balance governing the iml within a tissue, as described in supposition 1 
and constrained by postulate 1, obeys the integral equation 


3 tf P l dv = -[ p e v l -nda+ f (p e c - p\ + p^ - + p l s - p s t ) da. 

Ji 8 JSi 8 JB 

whose local form is the field equation 


(4.11) 


d\p l + p l divv l = ^ with = p l c - p c t + p* - p"l + pf - p\. 


(4.12) 


where is the density of mass supply to the interstitial miscible liquid, i.e., the 
rate per unit volume at which mass is being moved from the three other phases 
into the iml, minus the rate at which it is being removed from the iml and placed 
into the other three phases. From this equation, it immediately follows that 

trl/ = 0 (4.13) 


whenever the liquid phase is incompressible, which is a good assumption. 

Finally, the mass balance governing the pos within a tissue, as described in 
supposition 1 and constrained by postulate 1, obeys the integral equation 

9 1 f P s du = - f pV • n dr/ + f (p^ - Ps + Pi ~ Pi + Pm - PT) dv > ( 4 - 14 ) 
Jb Jsb Jb 

whose local form is the field equation 

d)p* + /div v s = c s with c s = pf + p* - 'p m s , (4.15) 

where c s is the density of mass supply to the population of solutes, i.e., the rate 
per unit volume at which mass is being moved from the three other phases into the 
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pos, minus the rate at which it is being removed from the pos and placed into the 
other three phases. From this equation, it immediately follows that 


tr V = 0 (4.16) 

whenever the solutes are incompressible, which is probably a reasonable assumption. 
A means by which the mass density for a pos can be quantified is via 

N s N s 

p s = J2 p s j = 2Z n i m j - ( 4 - 11 ) 

7 = 1 7=1 

where nj is the molar concentration (number of moles per cubic centimeter) and 
Mj is the molecular weight (grams per mole) of solute molecule j . 

The four mass supplies, which are mass rates per unit volume, satisfy the identity 

£>“ = 0, (4.18) 

in accordance with their definitions given in equations (4.5, 4.9, 4.12 & 4.15). Hence, 
our physical laws for tissue are a special case of Truesdell’s (1957) physical laws for 
chemically reacting mixtures with diffusion. 

From the identity in equation (4.3), along with equations (3.5, 3.7 & 3.14), 
the sum of the constituent mass-balance laws listed in equations (4.5, 4.9, 4.12 
& 4.15) equates with the conservation of mass stated in equation (4.2). In this 
regard, our theory differs from the existing theories for tissue growth discussed 
in the Introduction. Our theory, like Truesdell’s (1957) and Bowen’s (1969, 1976) 
theories, describes a closed-system continuum wherein mass is conserved both at a 
mass point and throughout the body. Prior tissue growth theories describe open- 
system continua wherein mass can be either created or destroyed at a mass point, 
with a tacit implication that mass is somehow being conserved at some higher level 
of organization such as the organism. 

It is convenient to express the balance equation that governs mass in each of 
the constituent phases as a single integral equation 

d t f p“dv = -f p a v a -nda + f <c“du, (4.19) 

is JSi 3 J£ 

whose local form is 

d“p“ + p“divi;“ = c a , (4.20) 

wherein a e {c, l , m, s}. In a general sense, the density of mass supply c a represents 
a sum of rates at which masses are leaving the three other phases to be absorbed 
by the a th phase, minus a sum of rates at which masses are departing from the a th 
phase to be absorbed by the other three phases. 

The equation for mass balance in the a th constituent, equation (4.20), and the 
associated formula governing mass balance in the overall mixture, equation (4.18), 
comprise what is theorem 1 in Truesdell’s (1957) monumental paper. 

A balance formula of the form 

d (V + p c div v m = div(£> grad p c ) + (4.21) 
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has been used by Tranquillo <fc Murray (1993) in their modeling of wound closure, 
and by Gomez-Benito et al. (2006) in their modeling of bone healing, with the 
mass supply <c c introducing terms akin to those of the logistics equation in order to 
model cell mitosis (birth) and apoptosis (death). Equation (4.21) is an alternative 
way to write equation (4.9) for mass balance in the poc, given the definition for v c 
provided in equation (3.13). 

(i) Rates for Mass- Averaged Fields 

Consider an arbitrary field /I defined as the mass-weighted average of its con- 
stituents 

P = or equivalently pfl — ^ p a . (4.22) 

From equation (3.18), it follows that 

d r/ 0“ = d“y6“ - grad /3“ • u a , (4.23) 

whilst from equations (3.15, 3.18 & 4.3), and an application of the chain rule to 
equation (3.7), viz., pd t p, a = d t p a — ix a d t p, which is further simplified via equation 
(4.2) for mass conservation, it follows that the balance law governing mass in each 
constituent phase, equation (4.20), can be recast as 

pd t fi a — — div(p“«“) + <c“, (4.24) 

such that equations (4.22-4.24), in the company of equation (4.3), collectively yield 
(cf. equation 1.2.17 in Bowen 1976) 

pd t (fi a p tt ) = p a d t p a + p tt pd tf i tt 

= p a d a t p a - div(p“ /?“«“) + 

and therefore 

pd t p = J2(p ad “P a -div(p a p a u a ) + 

which is a formula of extreme importance in mixture theory. This formula estab- 
lishes how local fields map into their associated global fields. It is the predominant 
averaging formula of mixture theory. 

As a matter of illustration, 

px = J](p a d>“ -div(p“i;“0 u f ) + (4.27) 

quantifies the bary centric acceleration of a mixture. 

5. Balance Laws for Tissues 

In the prior section, we have shown that tissues can be treated as a special case 
of chemically reacting mixtures with diffusible constituents. Consequently, one can 
immediately write down the physical laws that govern them. The balance laws for 
mass, momentum, and energy are from theorems 1-3 of Truesdell (1957). The field 
equations governing entropy are derived in the appendix. 

Tissues do not require the full generality of the theory for chemically react- 
ing mixtures with diffusion. Certain simplifying assumptions can be imposed; in 
particular: 


(4.25) 

(4.26) 
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Assumption 1. The temperatures of the separate phases are equal. 


Assumption 2. Second- and higher-order terms in diffusion velocity can be ne- 
glected. 

Assumption 3. All moments of momenta can be neglected. 

In accordance with these assumptions, the physical laws that govern the mass, 
momentum, energy, and entropy of each constituent in such a mixture are respec- 
tively 

d “p a = -p a div»“ + «“, 

p a d*v* =divT“ + p a g + p a , T“ = (T“) T , 
p a d“e“ - tr(T“D") - div?“ + p a r a + 
p a d?r?“ = -di y{q*/Q) + p a r a /9 + * a + y a , 

where a € {c,l,m,sj in our case. The laws that govern momentum and entropy 
utilized assumptions 1 & 3 in their derivations. Here p a , e“, r“, and rf are the mass, 
internal energy, heat production, and entropy densities of the a th constituent, while 
y a is its entropy production, 6 is temperature, g is the gravity vector (the sole body 
force considered here), v a and q a is the velocity and heat flux vectors of the a th 
constituent, and T“ is its partial stress. Quantities p a , e a . and are the rates 
at which momentum, internal energy, and entropy are being supplied to the a th 
phase from the other three phases through processes linked to some mass supply 
c a and/or some diffusion velocity 

Their counterparts, the physical laws that govern the overall mixture, are sat- 
isfied implicitly by the following set of constraint equations 


J2(p* + c«u«) 
+ tr (T“grad«“) + 

+ di v( P a f a u a /e) + ,cv) 


= 0 , 
= 0 . 
= 0 , 
= 0 . 
> 0 , 


(5.2) 


where assumptions 1 & 2 have been imposed, and where \[f a — e a — r] a 6 denotes 
the Helmholtz free energy of the a th phase. The first formula in equation (5.2) is 
just equation (4.18). The second formula, like the first, is as it appears in Truesdell 
(1957). However, the third formula is different, because of how work is considered 
to map (see equation 5.5). In our formulation, the additional work lost or gained 
through diffusion tr (T“grad;/“) is accounted for. In Truesdell’s (1957) and Bowen’s 
(1976) formulation, this effect is quantified by p a • which is associated with a 
different heat flux map than is given in equation (5.5). The remaining two formulae 
in equation (5.2) are new to the best of my knowledge, and are derived in the 
appendix. 
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The formulas in equation (5.2) are equivalent to the well-known field equations 


d t p — —p div x , 

pd t x — div T + pg, T = T t , 

pd t e = tr(TD) — div# + pr, 

pd t i ] — —d\\(q/Q) + pr/6 + y, y > 0, 


which are classic in their construction. In accordance with assumption 2, as it 
pertains to the internal energy, the mappings 


P = X>“’ * = £>“»“, e = ^ = (5-4) 

establish how the various local fields that reside on the left-hand sides of the for- 
mulae in equation (5.1), as arguments of their material derivatives, relate to their 
continuum counterparts in the formulae of equation (5.3). Also in accordance with 
assumption 2, which has been imposed on all of the formulae to follow except for 
that of the heat production, one obtains the mappings 


T = T“, 

tr(TD) = ^(tr(T“D“) — tr (T“grad */“)), 
q = ^(q 01 + p a e a u a ), 
r = Y,H a r a , 


(5-5) 


which establish how the various local fields relate to their global counterparts that 
reside on the right-hand sides of these formulae. 

See, for example, Bowen (1976) for the derivations and representations that hold 
whenever asumptions 1-3 do not apply. We point out that Bowen (1976) further 
decomposes tr(T“grad«“) into div(T“w“) — «“ ■ div T“, which is mathematically 
correct, and in doing so he is lead to different mapping formulae for q and e a 
than those that are given in equations (5.2 & 5.5). Both of our formulations are 
mathematically correct. We find our formulation to be more intuitive. Bowen (1976) 
refers to our q as k. 


6. Constitutive Structure 

The customary means by which constitutive equations are derived is to introduce 
a Lagendre transformation that swaps entropy for temperature as an independent 
thermodynamic variable, with temperature having the advantage of being capable 
of measurement. This particular transformation exchanges the internal energy with 
the Helmholtz free energy as the thermodynamic state function, which arose nat- 
urally in our derivation of entropy production, equation (A. 8), and allows one to 
combine the formulae for the first- and second-laws of thermodynamics into a single 
formula, with the outcome being a Gibbs-like equation. 
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Let the velocity fields for the four constituents of tissue be described by 
v e = 

V m = d t <p(X m ,t), 

N s 

grad to pi (6.1) 

i= 1 

N c , N s 

V c = v m - J2 ( 5)f grad In pf § rad ln 
i = l' 1 = 1 

in accordance with supposition 1. 

A common form adopted for the constitutive structure of momentum supply is 
(cf. equation 2.1.12 in Bowen 1976) 

p a = - <r“Vad p P -J2^ vt> ~ rgrad T, (6.2) 

P P 

where the first term on the right-hand side introduces a buoyancy effect, the second 
term introduces a Stokes drag effect, and the last term introduces a Soret thermal- 
diffusion effect, with <r“^, and £“ denoting material constants. For example, the 
biphasic theory of Klisch & Lotz (2000) utilizes Stokean drag as their constitutive 
equation for momentum supply in their modeling of the annulus fibrosus. 


Appendix. Balance of Entropy With Its Production 

Unlike axiom 1, which is a law of conservation (i.e. , an equality), the second law 
of thermodynamics in its classic presentation establishes an inequality. Following 
the approach promulgated by Truesdell & Noll (1965, pg. 295), one can postulate 
the existence of an entropy production term that ‘balances’ this classic entropy 
inequality, thereby putting the axiom for the second law into a format that can be 
readily applied down to the level of the sub-continua of a mixture. 

Axiom 2. The rate at which entropy increases inside of £ equals the flux of entropy 
entering across 8£, plus the entropy expended by the heat fluxing across its surface 
8d B and by the heat generating within its volume £, plus a non-negative rate of 
entropy generation created internally by other irreversible processes. 

Axiom 2 is a statement of the second law of rational thermodynamics, whose 
mathematical interpretation is 

3 f / p?)dv = — pr/x-nda— / {q/9)-nda+ / (pr/0 + y) du, (A.l) 

Jb Jsb Jsb Jb 
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which satisfies a Clausius-Duhem constraint via the integral inequality 

I y du > 0, (A. 2) 

Jb 

wherein rj and r are the entropy and rate of heat production per unit mass, 6 is the 
absolute temperature, q is the heat flux, and y is the rate of entropy production per 
unit volume brought about by internal irreversible processes. Following conventional 
arguments, and simplifying with equation (4.2) for mass conservation, these integral 
equations are equivalent to the field equations 

pd t q — -div(q/0) + pr/9 + y, (A. 3) 

and 

y > 0, (A. 4) 

where the latter formula handles the inequality from the second law of thermody- 
namics, as it applies to a continuum. These results are well known. 

Postulate 2. For a phase in a mixture, the rate at which entropy increases inside 
of 3 equals the flux of entropy entering across 83, plus the entropy expended 
by the heat fluxing across its surface 83 and by the heat generating within its 
volume 3, plus a non negative rate of entropy generation created internally by 
other irreversible processes, plus any entropies that are being exchanged between 
it and the other phases within 3. 

Given that equation (4.20) establishes the law governing mass balance within 
the phases of a mixture then, in accordance with postulate 2 and assumption 1, a 
balance of the entropy with its production in the a* 11 constituent requires that 

d t [ p a r] a dv — — [ p a if v a ■ n da — [ ( q a /9)-nda 
Jb Jsb Jsb 

+ J (p a r a /0 + y a + +«V) du, (A. 5) 

where rf is the entropy per unit mass and y a is the rate of entropy production per 
unit volume, both within the a th phase, while a“ is the density of entropy supply 
originating from an exchange of entropies between it and the other constituent 
phases. When written as a field equation, simplifying with equation (4.20) for mass 
balance, one gets 

p a d a t t] a = -6iv(q a /9) + p a r a /6 + y a + d“, (A.6) 

which is a balance law that governs entropy and its production at the constituent 
level. By introducing the notions of entropy production and entropy supply, we are 
able to construct an entropy law for the individual phases of a mixture in a manner 
that is consistent with the construction of their mass, momentum, and energy laws. 
In this regard, our approach differs from the classic approach outlined in Bowen 
(1976) wherein y, y a , and d“ have not been introduced. 

It is important to point out that postulate 2 does not carry over a ‘non-negative’ 
constraint on y a that is otherwise imposed on y in axiom 2; in other words, it is 
admissible for a y a to be negative as-long-as y is always non-negative. 
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By assigning the mappings 

y = J2y a and pq = ^p a q a , (A.7) 

one is able to express pd,q in terms of local fields via equation (4.26) leading to 

+ di v(p a xls a u a /0) + -cV) = 0, (A. 8) 

where the mappings for heat flux q and heat production r obtained from the en- 
ergy balance formulae have been applied, as listed in equation (5.5) wherein the 
expression for q requires an application of assumption 2. Here \j/ a — e a — q a 9 de- 
fines the Helmholtz free energy for the a th constituent. Equation (A. 8) forces the 
entropy supply to be in balance with the entropy lost or gained through diffusion 
and through mass supply. 

(a) Helmholtz Free Energy Formulation 

The conservation of energy listed in equation (5.3) and the balance of entropy 
with its production derived in equation (A. 3) combine to produce 

p d t x[r — —pi]d t 9 + tr(TD) — q ■ grad (In 9) — yd, (A. 9) 

where \fr — e — q9 is the Helmholtz free energy function. 

Likewise, the balance law for energy listed in equation (5.1) and the formula 
balancing entropy with its production derived in equation (A. 6), both of which 
pertain to the constituents of a mixture, and combine to produce 

p a d^xfr 01 = -A“ d “0 + tr(T“D“) - q a ■ grad (In 9) - y a 9 + a 01 , (A.10) 

where a a = e a — <x a 9 is the density of free- energy supply associated with the a th 
phase. Its constraint equation is obtained by multiplying temperature with the 
fourth formula in equation (5.2) and subtracting that from the third formula in 
equation (5.2). 
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